clear; close all; clc;

%% 固定忆阻器参数
a = 2;
b = 0.83;
c = 0.14;

%% 模拟总步数（建议足够长以达到稳态效果）
N = 1000;
n = 1:N;  % 离散时间索引

%% ========== 图(a) — 固定 ω=0.15, 初值 φ0=6, 不同 A = 0.8, 0.7, 0.6 ==========

omega_fixed = 0.15;
phi0_fixed = 6;
A_vals = [0.8, 0.7, 0.6];

figure;
subplot(1,3,1);  % 第1个子图
hold on;

for A = A_vals
    % 离散信号 v(n)
    v = A .* sin(omega_fixed .* n);

    % 预分配
    phi = zeros(1, N);
    i   = zeros(1, N);

    % 初值
    phi(1) = phi0_fixed;

    % 根据离散忆阻器模型迭代
    for k = 1:N-1
        i(k) = a * sin(b + c*phi(k)) * v(k);
        phi(k+1) = phi(k) + v(k);
    end
    i(N) = a * sin(b + c*phi(N)) * v(N);

    % 绘制电压v与电流i的关系曲线
    plot(i, v, 'LineWidth', 1.5);
end

title('(a): 固定 \omega=0.15, \phi_0=6, 不同A');
xlabel('输入电压 v');
ylabel('输出电流 i');
legend('A=0.8','A=0.7','A=0.6','Location','best');
grid on;
hold off;

%% ========== 图(b) — 固定 ω=0.15, A=0.8, 不同初值 φ0=6, 4, 2 ==========

A_fixed = 0.8;
phi0_vals = [6, 4, 2];

subplot(1,3,2);  % 第2个子图
hold on;

for phi0 = phi0_vals
    v = A_fixed .* sin(omega_fixed .* n);

    phi = zeros(1, N);
    i   = zeros(1, N);

    phi(1) = phi0;

    for k = 1:N-1
        i(k) = a * sin(b + c*phi(k)) * v(k);
        phi(k+1) = phi(k) + v(k);
    end
    i(N) = a * sin(b + c*phi(N)) * v(N);

    plot(i, v, 'LineWidth', 1.5);
end

title('(b): 固定 \omega=0.15, A=0.8, 不同 \phi_0');
xlabel('输入电压 v');
ylabel('输出电流 i');
legend('\phi_0=6','\phi_0=4','\phi_0=2','Location','best');
grid on;
hold off;

%% ========== 图(c) — 固定 A=0.8, φ0=6, 不同频率 ω=0.15, 0.2, 0.25 ==========

omega_vals = [0.15, 0.2, 0.25];
phi0_fixed = 6;
A_fixed = 0.8;

subplot(1,3,3);  % 第3个子图
hold on;

for omega = omega_vals
    v = A_fixed .* sin(omega .* n);

    phi = zeros(1, N);
    i   = zeros(1, N);

    phi(1) = phi0_fixed;

    for k = 1:N-1
        i(k) = a * sin(b + c*phi(k)) * v(k);
        phi(k+1) = phi(k) + v(k);
    end
    i(N) = a * sin(b + c*phi(N)) * v(N);

    plot(i, v, 'LineWidth', 1.5);
end

title('(c): 固定 A=0.8, \phi_0=6, 不同 \omega');
xlabel('输入电压 v');
ylabel('输出电流 i');
legend('\omega=0.15','\omega=0.2','\omega=0.25','Location','best');
grid on;
hold off;

%% 调整图形外观
set(gcf, 'Position', [100, 100, 1400, 400]);
